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ABSTRACT 

In this paper we extend the SIR epidemic model, considered in [1], from a deterministic framework to a stochastic 
one, and we formulate it as a stochastic differential equation. We present a Picard iteration method of proving verification 
theorem for the existence and uniqueness of the solution in global behavior. As a side effect of this study, we move one the 
development of the numerical treatment to solve the considered problem using the Milstein’s scheme. A comparative 
numerical study is done using some values of the intensity of the stochastic environment around the endemic equilibrium. 
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1. INTRODUCTION 

In 1927, 1932 and 1933 W.O. Kermack and A.G. McKendrick published a series of papers titled Contributions to 
the mathematical theory of epidemics. Those papers are often seen as the basis of further research in mathematical 
(especially deterministic) modelling of the spread of infectious diseases. The first three papers of Kermack and 
McKendrick are reprinted in [10]. In the book of Diekmann and Heesterbeek [11] some of the results of Kermack and 
McKendrick and other deterministic models are presented and explained. 

Other deterministic epidemiology models were then developed in papers by Ross, Ross and Hudson, Martini, and 
Lotka [5] [8] [9]. Mathematical epidemiology seems to have grown exponentially starting in the middle of the 20th century 
(the first edition in 1957 of Bailey’s book [7] is an important landmark), so that a tremendous variety of models have now 
been formulated, mathematically analyzed, and applied to infectious diseases. 

Allen [4] discusses a stochastic model of the above SIS epidemic model. However this is done by constructing a 
stochastic differential equation (SDE) approximation to the continuous time. The latter is obtained by assuming that events 
occurring at a constant rate in the deterministic model occur according to a Poisson process with the same rate. 
McCormack and Allen [12] construct a similar SDE approximation to an SIS multihost epidemic model and explore the 
stochastic and deterministic models numerically. Other previous work on parameter perturbation in epidemic models seems 
to have concentrated on the SIR model. Tornatore, Buccellato and Vetro [6] discuss an SDE SIR system with and without 
delay with a similar parameter perturbation. However this is not the only way to introduce stochasticity into the model. 

The population is divided into three distinct classes: the susceptible S, ’’healthy individuals who can catch the 
disease”, the infected I , ’’those who have the disease and can transmit it”, and the removed R, ’’individuals who have had 
the disease and are now immune to the infection” (or removed from further propagation of the dis- ease by some other 
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means). Schematically, the individual goes through consecutive states S — > I — > R. Such models are often called the SIR 
models. 


2. STOCHASTIC STUDY 


The SIR Model is used in the modeling of infectious diseases by computing the amount of people in a closed 
population that are susceptible, infected, or recovered at a given period of time. 


The model is also used by researchers and health officials to explain the increase and decrease in people needing 


medical care for a certain disease during an epidemic. The SIR model is the basis for other similar models. The SI model, 


once again. The common cold can be modeled with the SI model. There is also the SEIR model, where people are 
categorized as susceptible, exposed, infected, or recovered. The SIR model can be adjusted to include variation due to 
seasonal changes as seen by Bauch and Earn. [13] 

2.1. Stochastic Differential Equation for SIR Model 

we consider the SIR model given by : 


The model has a susceptible group designated by S, an infected group I, and a recovered group R with permanent 
immunity, rc is the intrinsic growth rate of susceptible, k is the carrying capacity of the susceptible in the absence of 
infective, a is the maximum values of per capita reduction rate of S du to I, a is half saturation constants, y is the natural 
recover rate from infection and 8 is the death rate of recovered populations. 

In the filtered probability space (Q, F , (Ft )t>0 , P), we consider the Stochastic version of deterministic SIR model 
(2.1), given by : 


with the initial conditions S(0) = SO , I (0) = 10 , R(0) = R0 , and S, I and R represent the number of susceptible, 
infective and recovered individuals, respectively. gS F (S, I , R), ol G(S, I , R) and oR H (S, I , R) are the diffusion 
coefficients with oi (i = S, R, I ) are real constants and known as the intensity of environmental fluctuations, Wi (i = S, R, 
I) are independent standard Brownian motions. 

The aim of the following subsection is to study the existence and uniqueness of considered problem (2.2), than to 


also known as the SIS model, is the model where once a person is no longer infectious, this per- son becomes susceptible 



( 2 . 1 ) 



( 2 . 2 ) 
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present the numerical resolution. Further- more we complete our study by giving numerical resolution for different values 
of ok with {k = S, I , R}. 

2.2. Existence and Uniqueness of Solution 
2.2.1. Global Behavior 

Throughout the rest of this chapter, let (Q, F, (Ft)t>0, P) be a complete filtered probability space. In this space, let 
f (.) and g(.) be bounded measurable functions mapping Rk into Rk and into the space of real k x n matrices, respectively. 
We now Consider the general n-dimensional stochastic system : 

dXt = f (t, Xt )dt + g(t, Xt )dWt (2.3) 

on t > 0 with initial value X (0) = X0 , the solution is denoted by X (t, X0 ) : Assume that f (t, 0) = g(t, 0) = 0 for 
all t > 0 Thus, we consider solutions to 


.A f =J£ 0 + f f(s,X a )ds+ ( 
jo Jo 


)dW , 


(2.4) 


By a solution to the SIDE described by (2.4) we mean a continuous. Ft -adapted process X (.) which satisfies (2.4) 
with probability one. 


If the following assumptions are satisfied 


• The functions /, g. are measurable. 

For every t 6 [O.T], x.x' € R", and there exists two constants KiJu such that, 

• ||/(i,*)— /(i, i'll + ||s(t,i) — ff(t, i'll! < JTi(|*-*'||) 

• ll/Mf + MM)II 2 < k?P + !M 2 } 

• Xo are square- integrable. 


Then the solution of the system (2.3) is unique, for more details see [3]. 


Remark 2.1. If we Taking X = (S.I,R), 


/ = Mi- 


q SI aSI 


k } 1 T al ’ 1 + al 


-iw 


g = I a s F (S, I, R),c,G(S, W), ”rH(S : I, R )\ , 
then the SDE (2.2) has an strong unique solution. 
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2-3. Proof 
Case 1: 

Let’s consider t £ [0, T ], 

X = {SJ,R), X' = (S\i\r‘) e Q =]si,s 2 [x]ii,i 2 [x]ri,r 3 [ with 

and 7*2 are the constants, and we defined |X|| = i/E St + P + R* 
We take F{S : L R) = S , G{S, /,i?) = / and (S, /, fl) = fl. 

We have g = (a $S, giI.orR). So: 


||/(f, X) - f(t, X')|| + ||g(t, X) - 5 (t. X')|| < K t (\\X - X'l) 
With 


| 9 (t,X)- 9 (i,J?')||<A- \\X -X' || 


With K = maxyogjO^a ^) | 

Let’s consider / = (/^ / 2 , / 3 ) With: 



1 + al 

h = l 1 




So 



For the second inequality we have: 
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ll/a(t,X)|l 2 < 


OLS2l2 
1 + Gl'l 


+ 7*2 


Let's consider 

A"' — max 

And 


{(r e s 2 ( 1 + j) +y 


fts 2 ? 2 \ / 0 * 2*2 


di i J V 1 H~ oil 


|[$(1. J 0II < °r) l X l 


+*)*} 


So we have 


l|/M)l | 2 + m,X)\f < K' + maxialclaim* 

< Ki+xim* 


Case 2: 


We take F(S. /, A} = exp (S), G(S, I. A) = exp(/}. and H(S, I. A ) = exp(A). 
We have g = (c s exp(S). a } exp(/), a R exp( A) j . 

We Know that: 


||rtt,x)-ffM')|| = 


\ 


3 2 5 

; ^ exp(S) - exp(S')j + cr| || exp(i) - exp{/')j + tr^exp(A) - exp(A') j 


using the following classical relation of exponential Function exp(y) > 1 4- y Vy E 
1. We get. 


|exp(5) -exp{S')|| < |]exp{S)|| 


1 - exp(S J - S ) 


exp(s^)||S - S'\\ 


by the same reasoning we obtain. 

||exp(/) -exp(J')|| < exp(i a )||/-/'|| 

|j exp(A) - exp{A')|| < exp(r 2 )||A - A'|| 

which shows that: 

ll/M) - /(i.Jf'HI + b(l, A') - g(t, X')|l < A'idl-V - x'l) 
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with 


*1 


m ax 


K (ri 


, 2j>s 2 . , 

_l 1- aj.a 


S-2 


(1 + aii ) 2 


)’(“’ Q (l + «i 1 ) 2 +7)’ K,! } 


such that ifu = ifina^^expt^), exp(s 2 ),exp(r 2 ) 

For the second inequality we have: 
l|j(*.*)ll a = o r |exp(25) + er 2 exp(2/) + o-^exp(2it) 


) 


< K 2 


exp(2S)^l + exp(— 2(S - /)) + exp(-2(S - R))^ 

!Xpf2S2) 1 + — 7 — - + — ^ 777 


< A' 2 exp(2s2) 

< i\ 2 exp(2s 2 ) 

< A" 2 exp(2s 2 ) 

< /f 2 exp(2s 2 ) 

< K?\\Xf 


(5-i + l) 2 (S-R+ 1) 


(S-/+l) 2 (5- A+ l) 2 (S-fl + 1) 2 + (S-/+1) 2 ' 
(5 - / + 1) 2 (5 - R + l) 2 + {S - I + 1) 2 (5 - R + l) 2 


(5- /+ 1) 2 (S- R+ 1) 2 + (S-R+ 1) 2 + (5- J + l) 2 


{s 2 — il + l)^( s 2 — r l T + 2 + 4 $2 + ||^"||' 


And we know that 

||/(f.X)|| 2 < K' 

< K[ 2 

So we obtain: 

+ || 9 (t,X)|l 2 < JC?(2+|SJC|| a ) 

So 

||/(f..Y)|| 2 + I|9( , W)I| 2 < *?(1 + ll*ll 2 ) 

With Aj 2 > 2A’J 2 

2.4. endemic equilibrium 

The deterministic SIR system admits a unique endemic equilibrium E” = 
(S’,r,ii*). with 

= 7(1 + Qf*) 
a 

akar c — 2r c a7 - a 2 A - + vA 
2 ;,. a 2 7 

r* = Zr 

5 
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A = [akar c — 2r c a 7 — a 2 k] 2 + 4 r r a 2 7 (fcar r — r c 7 ) 

For more details, refer to [14] and references therein. 

2.5. Numerical Simulation and Dynamics Comparison 

In this study, we will construct and implement numerical methods for solving “Stochastic Differential Equations” 

(SDEs) using MATLAB. An SDE is a differential equation in which one or more of the terms, and hence the 
solution itself, is a stochastic process. Studying stochastic processes requires a departure from the familiar deterministic 
setting of ordinary and partial differential equations, into a world where the evolution of a quantity has an inherent random 
component and where the expected behavior of this quantity can be described in terms of probability distributions. 

Depending on the mathematical model, we have a number of numerical techniques at our disposal. Using the 
Milstein’s higher order method in [2] to find the strong solution of system (2.2) with given initial value and the values of 
parameters. There can be many different approaches on deriving SDE model from the deterministic model, which may lead 
to different dynamical outcomes. In this sub- section, we compare the dynamics of stochastic differential equations with 
large diffusion terms. 

In the first time taking the coefficients F (S, I , R) = S, G(S, I , R) = I and H (S, I , R) = R, the corresponding 
discretization of model is : 

= Si + ^r c 3^1 - ^ Ji(el - 1)A* 

' Ii+l = L ' + “ 7/[ ) M + y *(1? - 

l ^+1 = ^*(i — tfJAt + 7 li&t — 

(2.5) 

where ^i , rji and ^i , i = 1, ..., n, are the Gaussian random variables N(0, 1). 

2.5.1. Numerical Results 

In this part, we give the obtained numerical results the Milstein’s higher order method implemented using Matlab 
software. The parameters used are represented in following Table. 

Table* I , Values of the parameters 


Parameter 

k r c a 

a 7 

s 

Value 

100 0.5 2.3 

1.49 0.25 

0.43 
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Figure 1: Solution of Stochastic Model 

Figure 1 shows that the susceptible S, infective I and removed R population will be oscillating by the noise but 
around E* = (1.0093, 2.1807, 1.2679) with the noise intensities respectively as = 0.12, ai = 0.1, ar = 0.21. 

For as = 2.72, ai = 1.89, ar = 1.51. 


Point d'equilibre (1 .0093,2.1 807,1 .2679) 



Figure 2: Solution of Stochastic Model 

From Figure 2, one can see that with increasing the noise intensities, the solutions of model (2.2) will be 
oscillating strongly around the endemic point 

E* = (1.0093, 2.1807, 1.2679) of model (2.6). 


Remark 2.2. The effects of the intensity of noise levels From figures 1 and 2, we can conclude that, when the 
intensity of noise is small, the stochastic model preserves the property of the global stability. In this case, we can ignore 
noise and use the deterministic model to approximate the population dynamics. However, the large intensity of noise can 
force the solutions of model (2.2) to oscillate strongly around the disease-free. In these cases, we cannot ignore the effect 
of noise, therefore, we cannot use deterministic model but stochastic model to describe the population dynamics. 

In this part we change the functions F, G and H in order to compare the fluctuation of System for the same values 
of ak , k = S, I , R. taking the coefficients F (S, I , R) = exp(S), G(S, I , R) = exp(I ) and H (S, I , R) = exp(R), the 
corresponding discretization of model is : 


Impact Factor (JCC): 2.0346 


NAAS Rating 3.19 


Numerical Solution of Stochastic Epidemiological Model: Case of Sir Model 


49 


S I+1 = 5, + |V C S, (i - f ) - At + exp(s.)'/*ki + Y - i)Af 

/,+! = / ( 4- - 7/, j Af + CTf t'xp(/j)>/Afiji + ^ exp(I t )(yj - l)Af 

L .%* = /?, f 1 - <5)At + Tr/i At + a R exp(^ ) + 4 ^p(^)(Cf - l)Af (2<6 - 

where , rji and , i = 1, n, are the Gaussian random variables N(0, 1). 

Taking the same values of as = 0.12, ai = 0.1, ar = 0.21 using in Figure 1, above in view to compare both figures. 


Point d'equilibre (1 .0093,2.1807,1 .2679) 



Figure 3: Solution of Stochastic Model 


Figure 3 illustrate fluctuations around the equilibrium point E* = (1.0093, 2.1807, 1.2679). In fact when varying 
functions F, G and H in order to enlarge the diffusion coefficients, then the fluctuation increases compare to figure 1. When 
we take the same values as = 2.72, ai = 1.89, ar = 1.51 as figure 2, we have: 


Point d'equilibre (1 .0093,2.1 807,1 .2679) 



In this case we remark that the fluctuations around the equilibrium point has increased greatly in this figure. 

3. CONCLUSIONS 


In this paper, we extended the SIR epidemic model considered in [1]. In particular we analyze the epidemic model 
in a stochastic environment. The interest of this study lies in two aspects. 

First, it presents the existence and uniqueness of solution using the Picard iteration method in the global behavior. 
Second, it performs some numerical simulations to illustrate the analytical results of stochastic model (2.2) by referring to 
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Milstein’s higher order method, and then a comparative numerical study is done for different value of a around the 
equilibrium epidemic. To study the effect of perturbation noise on the deterministic SIR model (2.1), we stochastically 
perturb model (2.6) with respect to white noise around its endemic equilibrium depending on the intensities of noise a. We 
conclude when the intensities of noise are not sufficiently large, the population of the stochastic model may be 
stochastically permanent around of equilibrium epidemic. 
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